##### This script runs the linear mixed-effects models reported in the Table OA3.

# The following set of models account for the possible BE and WE of LFPR only 
# (i.e. LFPR is the only contextul-level explanatory variable) and do account for compositional effects

# Model A: two levels with RE at the COUNTRY-ROUND level
Model_A_WE_BE_LFPR_only <- lmer(red_sup ~ employed + unemployed + hinc_2_cop + hinc_3_dif + hinc_4_vdif + female + age + isced2 + isced3 + isced4 + isced56 + stfeco + lrscale + 
                  mean_LFPR_y0 + dev_mean_LFPR_y0 + (1 | cntry_round), 
                  data = ess_4_ld_models, weights = pspwght)
summary(Model_A_WE_BE_LFPR_only)
check_collinearity(Model_A_WE_BE_LFPR_only)


# Model B: two levels with RE at the COUNTRY level
Model_B_WE_BE_LFPR_only <- lmer(red_sup ~ employed + unemployed + hinc_2_cop + hinc_3_dif + hinc_4_vdif + female + age + isced2 + isced3 + isced4 + isced56 + stfeco + lrscale + 
                  mean_LFPR_y0 + dev_mean_LFPR_y0 + (1 | cntry), 
                  data = ess_4_ld_models, weights = pspwght)
summary(Model_B_WE_BE_LFPR_only)
check_collinearity(Model_B_WE_BE_LFPR_only)


# Model C: three levels with RE at ROUND level (3) and then COUNTRY-ROUND level (2)
Model_C_WE_BE_LFPR_only <- lmer(red_sup ~ employed + unemployed + hinc_2_cop + hinc_3_dif + hinc_4_vdif + female + age + isced2 + isced3 + isced4 + isced56 + stfeco + lrscale + 
                  mean_LFPR_y0 + dev_mean_LFPR_y0 + (1 | cntry_round) + (1 | essround), 
                  data = ess_4_ld_models, weights = pspwght)
summary(Model_C_WE_BE_LFPR_only)
check_collinearity(Model_C_WE_BE_LFPR_only)


# Model D: three levels with RE at COUNTRY level (3) and then COUNTRY-ROUND level (2)
Model_D_WE_BE_LFPR_only <- lmer(red_sup ~ employed + unemployed + hinc_2_cop + hinc_3_dif + hinc_4_vdif + female + age + isced2 + isced3 + isced4 + isced56 + stfeco + lrscale + 
                  mean_LFPR_y0 + dev_mean_LFPR_y0 + (1 | cntry_round) + (1 | cntry), 
                  data = ess_4_ld_models, weights = pspwght)
summary(Model_D_WE_BE_LFPR_only)
check_collinearity(Model_D_WE_BE_LFPR_only)


# Model E: cross-classified model
Model_E_WE_BE_LFPR_only <- lmer(red_sup ~ employed + unemployed + hinc_2_cop + hinc_3_dif + hinc_4_vdif + female + age + isced2 + isced3 + isced4 + isced56 + stfeco + lrscale + 
                  mean_LFPR_y0 + dev_mean_LFPR_y0 + (1 | essround) + (1 | cntry), 
                  data = ess_4_ld_models, weights = pspwght)
summary(Model_E_WE_BE_LFPR_only)
check_collinearity(Model_E_WE_BE_LFPR_only)


# Model F: full model
Model_F_WE_BE_LFPR_only <- lmer(red_sup ~ employed + unemployed + hinc_2_cop + hinc_3_dif + hinc_4_vdif + female + age + isced2 + isced3 + isced4 + isced56 + stfeco + lrscale + 
                  mean_LFPR_y0 + dev_mean_LFPR_y0 + (1 | cntry_round) + (1 | essround) + (1 | cntry), 
                  data = ess_4_ld_models, weights = pspwght)
summary(Model_F_WE_BE_LFPR_only)
check_collinearity(Model_F_WE_BE_LFPR_only)


# Model G: two levels with RE at the YEAR level
Model_G_WE_BE_LFPR_only <- lmer(red_sup ~ employed + unemployed + hinc_2_cop + hinc_3_dif + hinc_4_vdif + female + age + isced2 + isced3 + isced4 + isced56 + stfeco + lrscale + 
                  mean_LFPR_y0 + dev_mean_LFPR_y0 + (1 | essround), 
                  data = ess_4_ld_models, weights = pspwght)
summary(Model_G_WE_BE_LFPR_only)
check_collinearity(Model_G_WE_BE_LFPR_only)


# Plotting the model results
tab_model(Model_A_WE_BE_LFPR_only, Model_B_WE_BE_LFPR_only, Model_C_WE_BE_LFPR_only, Model_D_WE_BE_LFPR_only, Model_E_WE_BE_LFPR_only, Model_F_WE_BE_LFPR_only, Model_G_WE_BE_LFPR_only, digits = 3, digits.p = 4, digits.re = 4)



#################################################################################
# Storing and exporting results of the seven models into a TABLE OA3 of the paper #
#################################################################################

# Creating an empty data.frame into which the results will be stored
Models_WE_BE_LFPR_only <- data.frame(c("Intercept", "", "employed", "", "unemployed", "", 
                               "hinc_2_cop", "", "hinc_3_dif", "", "hinc_4_vdif", "", 
                               "female", "", "age", "", "isced2", "", "isced3", "", 
                               "isced4", "", "isced56", "", "stfeco", "", "lrscale", "", 
                               "mean_LFPR_y0", "", "dev_mean_LFPR_y0", "", "Var_cntry", "Var_essround", 
                               "Var_cntry_round", "Var_individual", "ICC_cntry", "ICC_essround", "ICC_cntry_round"), matrix(data = rep(NA, 39*7), nrow = 39))
colnames(Models_WE_BE_LFPR_only) <- c("Parameter", "Model A", "Model B", "Model C", "Model D", "Model E", "Model F", "Model G")

# Using the programmed function to create Table OA3 of the article
Models_WE_BE_LFPR_only <- assign_fe_results(model_name = Model_A_WE_BE_LFPR_only, assignment_object = Models_WE_BE_LFPR_only, column_number = 2)
Models_WE_BE_LFPR_only <- assign_fe_results(model_name = Model_B_WE_BE_LFPR_only, assignment_object = Models_WE_BE_LFPR_only, column_number = 3)
Models_WE_BE_LFPR_only <- assign_fe_results(model_name = Model_C_WE_BE_LFPR_only, assignment_object = Models_WE_BE_LFPR_only, column_number = 4)
Models_WE_BE_LFPR_only <- assign_fe_results(model_name = Model_D_WE_BE_LFPR_only, assignment_object = Models_WE_BE_LFPR_only, column_number = 5)
Models_WE_BE_LFPR_only <- assign_fe_results(model_name = Model_E_WE_BE_LFPR_only, assignment_object = Models_WE_BE_LFPR_only, column_number = 6)
Models_WE_BE_LFPR_only <- assign_fe_results(model_name = Model_F_WE_BE_LFPR_only, assignment_object = Models_WE_BE_LFPR_only, column_number = 7)
Models_WE_BE_LFPR_only <- assign_fe_results(model_name = Model_G_WE_BE_LFPR_only, assignment_object = Models_WE_BE_LFPR_only, column_number = 8)
View(Models_WE_BE_LFPR_only)

### Manual assignment of variance components and ICCs
# Model A: RE at cntry_round
Models_WE_BE_LFPR_only[Models_WE_BE_LFPR_only$Parameter == "Var_cntry_round", "Model A"] <- as.data.frame(VarCorr(Model_A_WE_BE_LFPR_only))[, "vcov"][1]
Models_WE_BE_LFPR_only[Models_WE_BE_LFPR_only$Parameter == "Var_individual", "Model A"] <- as.data.frame(VarCorr(Model_A_WE_BE_LFPR_only))[, "vcov"][2]
icc(Model_A_WE_BE_LFPR_only, by_group = TRUE)
Models_WE_BE_LFPR_only[Models_WE_BE_LFPR_only$Parameter == "ICC_cntry_round", "Model A"] <- (as.data.frame(VarCorr(Model_A_WE_BE_LFPR_only))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_A_WE_BE_LFPR_only))[, "vcov"]))*100

# Model B: RE at cntry
Models_WE_BE_LFPR_only[Models_WE_BE_LFPR_only$Parameter == "Var_cntry", "Model B"] <- as.data.frame(VarCorr(Model_B_WE_BE_LFPR_only))[, "vcov"][1]
Models_WE_BE_LFPR_only[Models_WE_BE_LFPR_only$Parameter == "Var_individual", "Model B"] <- as.data.frame(VarCorr(Model_B_WE_BE_LFPR_only))[, "vcov"][2]
icc(Model_B_WE_BE_LFPR_only, by_group = TRUE)
Models_WE_BE_LFPR_only[Models_WE_BE_LFPR_only$Parameter == "ICC_cntry", "Model B"] <- (as.data.frame(VarCorr(Model_B_WE_BE_LFPR_only))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_B_WE_BE_LFPR_only))[, "vcov"]))*100

# Model C: RE at essround and cntry_round
Models_WE_BE_LFPR_only[Models_WE_BE_LFPR_only$Parameter == "Var_essround", "Model C"] <- as.data.frame(VarCorr(Model_C_WE_BE_LFPR_only))[, "vcov"][2]
Models_WE_BE_LFPR_only[Models_WE_BE_LFPR_only$Parameter == "Var_cntry_round", "Model C"] <- as.data.frame(VarCorr(Model_C_WE_BE_LFPR_only))[, "vcov"][1]
Models_WE_BE_LFPR_only[Models_WE_BE_LFPR_only$Parameter == "Var_individual", "Model C"] <- as.data.frame(VarCorr(Model_C_WE_BE_LFPR_only))[, "vcov"][3]
icc(Model_C_WE_BE_LFPR_only, by_group = TRUE)
Models_WE_BE_LFPR_only[Models_WE_BE_LFPR_only$Parameter == "ICC_essround", "Model C"] <- (as.data.frame(VarCorr(Model_C_WE_BE_LFPR_only))[, "vcov"][2]/sum(as.data.frame(VarCorr(Model_C_WE_BE_LFPR_only))[, "vcov"]))*100
Models_WE_BE_LFPR_only[Models_WE_BE_LFPR_only$Parameter == "ICC_cntry_round", "Model C"] <- (as.data.frame(VarCorr(Model_C_WE_BE_LFPR_only))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_C_WE_BE_LFPR_only))[, "vcov"]))*100

# Model D: three levels with RE at COUNTRY level (3) and then COUNTRY-ROUND level (2)
Models_WE_BE_LFPR_only[Models_WE_BE_LFPR_only$Parameter == "Var_cntry", "Model D"] <- as.data.frame(VarCorr(Model_D_WE_BE_LFPR_only))[, "vcov"][2]
Models_WE_BE_LFPR_only[Models_WE_BE_LFPR_only$Parameter == "Var_cntry_round", "Model D"] <- as.data.frame(VarCorr(Model_D_WE_BE_LFPR_only))[, "vcov"][1]
Models_WE_BE_LFPR_only[Models_WE_BE_LFPR_only$Parameter == "Var_individual", "Model D"] <- as.data.frame(VarCorr(Model_D_WE_BE_LFPR_only))[, "vcov"][3]
icc(Model_D_WE_BE_LFPR_only, by_group = TRUE)
Models_WE_BE_LFPR_only[Models_WE_BE_LFPR_only$Parameter == "ICC_cntry", "Model D"] <- (as.data.frame(VarCorr(Model_D_WE_BE_LFPR_only))[, "vcov"][2]/sum(as.data.frame(VarCorr(Model_D_WE_BE_LFPR_only))[, "vcov"]))*100
Models_WE_BE_LFPR_only[Models_WE_BE_LFPR_only$Parameter == "ICC_cntry_round", "Model D"] <- (as.data.frame(VarCorr(Model_D_WE_BE_LFPR_only))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_D_WE_BE_LFPR_only))[, "vcov"]))*100

# Model E: RE at country and ess_round
Models_WE_BE_LFPR_only[Models_WE_BE_LFPR_only$Parameter == "Var_cntry", "Model E"] <- as.data.frame(VarCorr(Model_E_WE_BE_LFPR_only))[, "vcov"][1]
Models_WE_BE_LFPR_only[Models_WE_BE_LFPR_only$Parameter == "Var_essround", "Model E"] <- as.data.frame(VarCorr(Model_E_WE_BE_LFPR_only))[, "vcov"][2]
Models_WE_BE_LFPR_only[Models_WE_BE_LFPR_only$Parameter == "Var_individual", "Model E"] <- as.data.frame(VarCorr(Model_E_WE_BE_LFPR_only))[, "vcov"][3]
icc(Model_E_WE_BE_LFPR_only, by_group = TRUE)
Models_WE_BE_LFPR_only[Models_WE_BE_LFPR_only$Parameter == "ICC_cntry", "Model E"] <- (as.data.frame(VarCorr(Model_E_WE_BE_LFPR_only))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_E_WE_BE_LFPR_only))[, "vcov"]))*100
Models_WE_BE_LFPR_only[Models_WE_BE_LFPR_only$Parameter == "ICC_essround", "Model E"] <- (as.data.frame(VarCorr(Model_E_WE_BE_LFPR_only))[, "vcov"][2]/sum(as.data.frame(VarCorr(Model_E_WE_BE_LFPR_only))[, "vcov"]))*100

# Model F: RE at all three effects
Models_WE_BE_LFPR_only[Models_WE_BE_LFPR_only$Parameter == "Var_cntry", "Model F"] <- as.data.frame(VarCorr(Model_F_WE_BE_LFPR_only))[, "vcov"][2]
Models_WE_BE_LFPR_only[Models_WE_BE_LFPR_only$Parameter == "Var_essround", "Model F"] <- as.data.frame(VarCorr(Model_F_WE_BE_LFPR_only))[, "vcov"][3]
Models_WE_BE_LFPR_only[Models_WE_BE_LFPR_only$Parameter == "Var_cntry_round", "Model F"] <- as.data.frame(VarCorr(Model_F_WE_BE_LFPR_only))[, "vcov"][1]
Models_WE_BE_LFPR_only[Models_WE_BE_LFPR_only$Parameter == "Var_individual", "Model F"] <- as.data.frame(VarCorr(Model_F_WE_BE_LFPR_only))[, "vcov"][4]
icc(Model_F_WE_BE_LFPR_only, by_group = TRUE)
Models_WE_BE_LFPR_only[Models_WE_BE_LFPR_only$Parameter == "ICC_cntry", "Model F"] <- (as.data.frame(VarCorr(Model_F_WE_BE_LFPR_only))[, "vcov"][2]/sum(as.data.frame(VarCorr(Model_F_WE_BE_LFPR_only))[, "vcov"]))*100
Models_WE_BE_LFPR_only[Models_WE_BE_LFPR_only$Parameter == "ICC_essround", "Model F"] <- (as.data.frame(VarCorr(Model_F_WE_BE_LFPR_only))[, "vcov"][3]/sum(as.data.frame(VarCorr(Model_F_WE_BE_LFPR_only))[, "vcov"]))*100
Models_WE_BE_LFPR_only[Models_WE_BE_LFPR_only$Parameter == "ICC_cntry_round", "Model F"] <- (as.data.frame(VarCorr(Model_F_WE_BE_LFPR_only))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_F_WE_BE_LFPR_only))[, "vcov"]))*100

# Model G: RE at essround
Models_WE_BE_LFPR_only[Models_WE_BE_LFPR_only$Parameter == "Var_essround", "Model G"] <- as.data.frame(VarCorr(Model_G_WE_BE_LFPR_only))[, "vcov"][1]
Models_WE_BE_LFPR_only[Models_WE_BE_LFPR_only$Parameter == "Var_individual", "Model G"] <- as.data.frame(VarCorr(Model_G_WE_BE_LFPR_only))[, "vcov"][2]
icc(Model_G_WE_BE_LFPR_only, by_group = TRUE)
Models_WE_BE_LFPR_only[Models_WE_BE_LFPR_only$Parameter == "ICC_essround", "Model G"] <- (as.data.frame(VarCorr(Model_G_WE_BE_LFPR_only))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_G_WE_BE_LFPR_only))[, "vcov"]))*100

# Exporting Table OA3 to a csv file
write.csv(x = Models_WE_BE_LFPR_only, file = "table_OA3_Models_BE_WE_LFPR_only.csv")



### Collecting the coefficient estimates for the LFPR BE contextual variable and storing them to a data.frame
only_LFPR_BE <- data.frame(c("Model A", "Model B", "Model C", "Model D", "Model E", "Model F", "Model G"),
                      matrix(data = rep(NA, 7*3), nrow = 7))
colnames(only_LFPR_BE) <- c("Model type", "Estimate", "Std. Error", "df")

only_LFPR_BE[only_LFPR_BE$`Model type` == "Model A", 2:4] <- coef(summary(Model_A_WE_BE_LFPR_only))[15, c("Estimate", "Std. Error", "df")]
only_LFPR_BE[only_LFPR_BE$`Model type` == "Model B", 2:4] <- coef(summary(Model_B_WE_BE_LFPR_only))[15, c("Estimate", "Std. Error", "df")]
only_LFPR_BE[only_LFPR_BE$`Model type` == "Model C", 2:4] <- coef(summary(Model_C_WE_BE_LFPR_only))[15, c("Estimate", "Std. Error", "df")]
only_LFPR_BE[only_LFPR_BE$`Model type` == "Model D", 2:4] <- coef(summary(Model_D_WE_BE_LFPR_only))[15, c("Estimate", "Std. Error", "df")]
only_LFPR_BE[only_LFPR_BE$`Model type` == "Model E", 2:4] <- coef(summary(Model_E_WE_BE_LFPR_only))[15, c("Estimate", "Std. Error", "df")]
only_LFPR_BE[only_LFPR_BE$`Model type` == "Model F", 2:4] <- coef(summary(Model_F_WE_BE_LFPR_only))[15, c("Estimate", "Std. Error", "df")]
only_LFPR_BE[only_LFPR_BE$`Model type` == "Model G", 2:4] <- coef(summary(Model_G_WE_BE_LFPR_only))[15, c("Estimate", "Std. Error", "df")]

# Computing the 95% confidence interval (based on student?s t-distribution)
only_LFPR_BE$lower_bound <- only_LFPR_BE$Estimate - only_LFPR_BE$`Std. Error`*qt(p = 0.975, df = only_LFPR_BE$df)
only_LFPR_BE$upper_bound <- only_LFPR_BE$Estimate + only_LFPR_BE$`Std. Error`*qt(p = 0.975, df = only_LFPR_BE$df)


### Collecting the coefficient estimates for the LFPR WE contextual variable and storing them to a data.frame
only_LFPR_WE <- data.frame(c("Model A", "Model B", "Model C", "Model D", "Model E", "Model F", "Model G"),
                      matrix(data = rep(NA, 7*3), nrow = 7))
colnames(only_LFPR_WE) <- c("Model type", "Estimate", "Std. Error", "df")

only_LFPR_WE[only_LFPR_WE$`Model type` == "Model A", 2:4] <- coef(summary(Model_A_WE_BE_LFPR_only))[16, c("Estimate", "Std. Error", "df")]
only_LFPR_WE[only_LFPR_WE$`Model type` == "Model B", 2:4] <- coef(summary(Model_B_WE_BE_LFPR_only))[16, c("Estimate", "Std. Error", "df")]
only_LFPR_WE[only_LFPR_WE$`Model type` == "Model C", 2:4] <- coef(summary(Model_C_WE_BE_LFPR_only))[16, c("Estimate", "Std. Error", "df")]
only_LFPR_WE[only_LFPR_WE$`Model type` == "Model D", 2:4] <- coef(summary(Model_D_WE_BE_LFPR_only))[16, c("Estimate", "Std. Error", "df")]
only_LFPR_WE[only_LFPR_WE$`Model type` == "Model E", 2:4] <- coef(summary(Model_E_WE_BE_LFPR_only))[16, c("Estimate", "Std. Error", "df")]
only_LFPR_WE[only_LFPR_WE$`Model type` == "Model F", 2:4] <- coef(summary(Model_F_WE_BE_LFPR_only))[16, c("Estimate", "Std. Error", "df")]
only_LFPR_WE[only_LFPR_WE$`Model type` == "Model G", 2:4] <- coef(summary(Model_G_WE_BE_LFPR_only))[16, c("Estimate", "Std. Error", "df")]

# Computing the 95% confidence interval (based on student?s t-distribution)
only_LFPR_WE$lower_bound <- only_LFPR_WE$Estimate - only_LFPR_WE$`Std. Error`*qt(p = 0.975, df = only_LFPR_WE$df)
only_LFPR_WE$upper_bound <- only_LFPR_WE$Estimate + only_LFPR_WE$`Std. Error`*qt(p = 0.975, df = only_LFPR_WE$df)



# Deleting model results from the workspace (as they are no longer needed)
rm(Model_A_WE_BE_LFPR_only, Model_B_WE_BE_LFPR_only, Model_C_WE_BE_LFPR_only, Model_D_WE_BE_LFPR_only, Model_E_WE_BE_LFPR_only, Model_F_WE_BE_LFPR_only, Model_G_WE_BE_LFPR_only)
